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We study the mean-field phase diagram of the repulsive Hubbard model in the Lieb lattice. 
Far from half-filling, the most stable phases are paramagnetism for low on-site interaction U/t 
and ferromagnetism for high U/t, as in the case of the mean-field phase diagram of the square 
lattice Hubbard model obtained by Dzierzawa [T. At half-filling, the ground state was found 
to be ferrimagnetic [a (n,n) spiral phase], in agreement with a theorem by Lieb [ 2 ]. The total 
magnetization approaches Lieb’s prediction as U/t becomes large. As we move away from half¬ 
filling, this ferrimagnetic phase becomes a (qi,qi) spiral phase with 51 m tv and then undergoes a 
series of first-order phase transitions, (51,51) —► (51,52) —» (51,0), with 52 « n/2, before becoming 
ferromagnetic at large U/t or paramagnetic at low U/t. 


I. INTRODUCTION 

The Hubbard model is one of the most studied models 
in the area of strongly correlated electron systems m- 
However, it remains unsolved for dimensionality larger 
than one. For the one-dimensional (ID) case, the exact 
solution is given by the Bethe Ansatz [Sj, while in the 
case of two dimensions (2D), the solution is known only 
in some limiting cases or by means of approximations, 
such as mean-field. The fermionic Hubbard model in a 
square lattice has long been known to display antiferro¬ 
magnetism (AF) at half-filling [5;. However, away from 
half-filling, the ground state magnetic ordering is still an 
open problem [7]- 

Extensions of the Hubbard model to 2D decorated lat¬ 
tices also show interesting features, such as flat band fer¬ 
romagnetism (F) [2] and Dirac cones [Sj. These decorated 
2D lattices fall into three classes: Lieb’s 0, Mielke’s 
[9] and Tasaki’s |TD]. The pursuit for metallic ferro¬ 
magnetism has motivated the search of crystal structures 
matching these decorated lattices. However, there are ex¬ 
perimental obstacles, such as the lifting of the flat-band 
degeneracy by the Jahn-Teller effect or the difficulty in 
controlling the filling of the lattice. An alternative exper¬ 
imental approach is to study quantum dot arrays with 
these geometries m- Decorated lattices can also be re¬ 
alized by manipulating cold atoms in optical lattices |12| . 

Here, we study one example of a 2D decorated lattice, 
the Lieb lattice, i.e., a line-centered square lattice m 
This kind of lattice can be obtained from the usual 2D 
square lattice by removing a quarter of its atoms (see 
Fig. 0- Each unit cell contains one atom of each kind: 
A, B and C. As a matter of fact, real materials can have 
their atoms arranged in a fashion resembling the Lieb 
lattice. Examples include the well-known high-T c super¬ 
conductors with weakly coupled Cu02 planes 0 0, 
such as La2- x Sr :E Cu04 and YBa2Cu30y, which can be 
studied using the perovskite lattice, a three-dimensional 
(3D) generalization of the Lieb lattice [IB] . 

Exact results for magnetism in the Lieb lattice are 
known. For example, an important theorem proven by 
Lieb [2] states that bipartite lattices (lattices with two 
sublattices, A and B, such that each site on sublattice A 
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FIG. 1 . A square lattice can be divided into four sublattices 
A, B, C and D. The circles represent atomic nuclei and the 
arrows represent spins. The Lieb lattice can be obtained by 
removing one of the sublattices. 


has its nearest neighbors on sublattice B, and vice versa) 
whose unit cell contains a different number of each kind 
of atom, have ferromagnetic ground states at half-filling. 
This is the case of the Lieb lattice BZI. as each unit cell 
contains one A atom and two B-like atoms. One com¬ 
mon argument is that these states are in fact ferrimag¬ 
netic )18| , in the sense that although each sublattice is 
ferromagnetic, the full lattice is antiferromagnetic, but 
the magnetization is finite due to the different number of 
atoms in each sublattice. This contrasts with the antifer¬ 
romagnetic ordering of the square lattice Hubbard model 
in this limit. Note that Lieb’s theorem only mentions the 
total magnetization per unit cell, not on-site magnetiza¬ 
tion amplitudes, which can be calculated using numerical 
methods, such as mean-field. This has been done for the 
multi-layer Lieb optical lattice at half-filling ma¬ 
in this work, we use a mean-field approach to com¬ 
pute the magnetic phase diagram of the Lieb lattice as a 
function of the average electron density n and Hubbard 
interaction U, thus going away from both half-filling and 
the tight-binding limit. The allowed magnetic phases are 
paramagnetism and spin spiral phases [20] . Ferro- and 
ferrimagnetism can be obtained as particular cases of spi¬ 
ral phases. Note that we do not consider spatial phase 
separation. In order to find such regions in the phase 
diagram, one needs to use the chemical potential as an 
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FIG. 2. (a) Plot of the tight-binding dispersion relation of 
the Hubbard model in the Lieb lattice and (b) the respective 
particle density of each sublattice, A, B or C, as a function of 
the total particle density. 


independent variable rather than using the par¬ 

ticle density. 

The tight-binding Hamiltonian of the Lieb lattice, H tl 
is given by [ 23 ] 

t E E [(4 ,yB x ,y + A.. y C,. V + H.C.) 

X=1 y=l K 1 / 

+ + 4 , y Cx-l,y + H.C.)] . 

L x ( Ly ) is the number of unit cells along the x (y) direc¬ 
tion. The hopping terms in the first line are intra-unit 
cell and the remaining ones are inter-unit cell. The eigen¬ 
values of H t originate three energy bands, one of which is 
flat. The dispersion relation for periodic boundary con¬ 
ditions is 


£± = ±2t^Jcos 2 y +COS 2 y, (2) 

for the non-flat energy bands, where k a = 27 Tn a /L a with 
n a = 0,1, • • ■ , L a and a £ {x, y}. The flat band is L x L y - 
fold degenerate with zero energy. The one-particle local¬ 
ized states associated with the flat bands can be written 
as 


shares, while the lower and upper bands involve all three 
sublattices A, B and C. 

The particle density of a sublattice equals the number 
of electrons at that sublattice divided by the number of 
atoms the sublattice comprises, or the number of unit 
cells, 


n * = zX 
n s = zrfc 


(4) 


nc 


N c 
L X Lnj 


In the non-interacting limit and at zero temperature, the 
particle density on sublattices A, B and C as a function 
of the global particle density (number of electrons in the 
whole lattice divided by the total number of sites) is as 
plotted in Fig. [2b] The plot can be interpreted as follows. 
Half of the probability density of the states in the lower 
dispersive band correspond to the sublattice A, while the 
other half is evenly distributed among sublattices B and 
C. Therefore, starting at n = 0, as we insert electrons in 
the system, half “choose” sublattice A, while the other 
half go to sublattices B or C. At n = 2/3, we reach the 
flat band at e = 0. At this point, all sites A are singly 
occupied, while sites B and C are quarter-filled. Any 
newly-added electrons will only go to sublattices B or C, 
because the flat band only comprises these two kinds of 
atoms and going to sublattice A would imply going to 
the upper dispersive band, which would lead to higher 
total energy. At n = 4/3, the flat band is completely 
filled, so that for n > 4/3 electrons occupy the upper 
dispersive band going to sites A or B/C at a ratio of 
2:1, as in the lower dispersive band, up to the maximum 
filling tia = tib = nc = 2. 

In this work, we address magnetism in the Lieb lat¬ 
tice by considering a finite on-site Coulomb repulsion U 
using a mean-field approach, and build a n — U phase di¬ 
agram. In the case of a square lattice, one assumes that 
the occupation number is the same in the whole lattice. 
Here, in the case of the Lieb lattice, we assume that the 
occupation number on each sublattice is the same as in 
the tight-binding limit, for any U (see Fig. 2b). This is 
the correct assumption for small U/t. Moreover, for large 
U/t , the results of Fig. [5] remain qualitatively the same 
for tia = tib = nc = n. 


Iloc; x,y) = ^ (Bl y - C^ y + 4, y -i - c l-i, y ) l vac > • 

(3) 

These states form a non-orthogonal basis of the flat band 
subspace. 

The three tight-binding energy bands of the Lieb lat¬ 
tice energy bands are shown in Fig. |2a| At the point 
(k x ,k v ) = (7r,7r), the three branches touch each other. 
Expanding the dispersion relation in Eq. [2] around this 
point, we find the Dirac cones e 2 = t 2 (fc 2 + k y ). The flat 
band is built up from B- and C-type orbitals in equal 


II. CALCULATIONS 

The interaction term of the Hubbard Hamiltonian is 

Hu = uY: n T n U ( 5 ) 

sites 

that is, the on-site Coulomb repulsion U times the num¬ 
ber of double occupancies in the lattice. Applying the 
mean-field approximation to the Hubbard Hamiltonian 
gives single-particle energies given by the eigenvalues of 
the 6x6 single-particle Hamiltonian Hmf [2 [mug, 
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plus the diagonal term 



/o 2 2 2 2 \ 

(ora — n A — n B — n c ). 


of U (the on-site interaction), we can control the ratio 
U/t , for example by applying pressure on the sample. 


This is a generalization of the Hamiltonian obtained in 
previous studies of the 2D square lattice Hubbard model 
(II HZ5112B] , which did not allow for different occupations 
in the sublattices. The magnetic phase of the system is 
defined by two order parameters: the vector q and the 
number to, as in the works by Dzierzawa [T] and Singh 
I25| . The vector q = {q x -,q y ) defines the orientation of 
the spins. For example, q x = 0 is a ferromagnetic phase 
along the x direction, q y = tt represents antiferromag¬ 
netism along the y direction, and other values of q x or 
q y give spin spiral phases. The paramagnetic phase is 
^-degenerate and is characterized by zero magnetization 
amplitude. The magnetization amplitude to can be iden¬ 
tified as the amplitude of the spin spiral wave, 

(Sf) = (cos(g- r),sin(g- f),0), (8) 

and appears during the mean-field calculations, when 
computing averages such as 

MW = (Si ) = (Sa,x + iS A ,y) = y , (9) 

for sublattice A. Fig. [3] shows what the configuration of 
the Lieb lattice looks like when q = (7r, 7r) and m is finite. 
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FIG. 3. According to our definition of q (the change in spin 
orientation between two consecutive lattice sites), a Lieb lat¬ 
tice with q = (n, it) is ferromagnetic within each sublattice. 
Moreover, because the magnetization amplitude m is the same 
in every site and each unit cell has two spins in the same di¬ 
rection and one spin in the opposite direction, the total spin 
per unit cell is nonzero, as predicted by Lieb [2i. 


From this point forward, we consider t = 1, so that U is 
given in units of t. It is important to remark that, exper¬ 
imentally, although we cannot directly control the value 


III. RESULTS AND DISCUSSION 

The n — U phase diagram is computed in the following 
way. For each point ( n,U ), the number of electrons N 
is well defined, so that we can add the lowest N mean- 
field energies and find the total energy of the system. 
By numerically minimizing this total energy (using, for 
instance, the algorithm in Ref. m) with respect to q x , 
q y and to, we find the values of these three magnetic 
order parameters which lead to the ground state for this 
pair (n, U). Repeating this process for all desired pairs, 
one obtains the phase diagram. 


A. Magnetization for high U/t 

The plot in Fig. [4] shows the mean-field ground state 
magnetization amplitude m as a function of n and U. 
This result is similar to that of the square lattice in most 
regions of the diagram. Indeed, for high U, the magneti¬ 
zation is proportional to n between n = 0 and n = 1, and 
proportional to 2 — n between n = 1 and n = 2, reflecting 
particle-hole symmetry. 

This proportionality can be justified analytically in the 
following way. For very high U, the tight-binding terms 
of the mean-field Hamiltonian given by Eq. [6] are but 
a small perturbation, which can be neglected as a first 
approximation. In this case, all sublattices become equiv¬ 
alent, implying n A = ns = «c = n, with t = 0. Diag¬ 
onalizing the new Hamiltonian gives two flat bands. A 
three-fold degenerate band at (f(n — to) and a three-fold 
degenerate band at ^ (n + to) . Distributing the electrons 
in the bands and adding up their energies, one obtains 
the total energy of the system by adding the diagonal 
term 3U ^ uc (jn 2 — n 2 ) (see Eq. [7]), so that having positive 
to yields the same bands as negative m. Let us assume 
in > 0. For n £ [0,1], electrons occupy only the lowest 
(degenerate) energy band, with energy ^(n — to). The 
total energy is then given by 

3U ^ UC ( to 2 - n 2 ) + ^ J2( n ~ m )- ( 10 ) 

N 

The result of the summation is simply N(n — to) = 
3 nL uc (n — to). Minimizing this with respect to to gives 
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FIG. 4. Mean-field ground state magnetization (m) of the 
Lieb lattice as a function of n and U. The plot is very similar 
to that of the 2D square lattice. The most noticeable differ¬ 
ence between the two is that, while the square lattice has zero 
m in the vicinity of the point ( n , U) = (1, 0), the Lieb lattice 
has finite m in this region of the diagram, more specifically 
between n = 2/3 and n = 4/3. 


the expected result to = n. Performing an analogous cal¬ 
culation assuming n > 1 yields the relation m = 2 — n, 
i.e., the other half of the plot in Fig. [4] for high U. 

B. Magnetization for low U/t 


not depend on to, and having finite to would only increase 
the total energy due to the diagonal term in Eq. [7] As 
the total particle density reaches n ~ 2/3 (getting closer 
to 2/3 as U approaches 0), we start to fill the nearly flat 
bands at e = 0. This is the point at which a finite to can 
be used to lower the energy of one of the flat bands, thus 
lowering the total energy of the system. After the lower 
flat band has been filled (note that finite U induces some 
modulation of the flat bands, but the argument is valid 
for small perturbations), we are at n = 1 and start filling 
the upper flat band. Now, it becomes advantageous to 
lower the value of to, so as to reduce the energy of this 
band. Finally, at n k 4/3, only the two higher-energy 
dispersive bands remain empty and between n = 4/3 
and n = 2, the value of to goes back to zero, for the same 
reason as when filling the two lowest-energy bands. 

Let us now compare these assertions with our numeri¬ 
cal results in Fig. [4] At small U and far from half-filling 
(outside the n interval [2/3;4/3]), the ground state of 
the system is paramagnetic (to = 0 ), coinciding with the 
square lattice result (see Fig. [5a]). On the other hand, 
inside the interval n £ [2/3; 4/3J7 the square lattice be¬ 
comes paramagnetic (except at exactly n = 1, where it 
is antiferromagnetic, and in a very small region around 
n=l, where a spiral phase arises; the width of this re¬ 
gion goes to zero as U/t —> oo) while our result suggests 
that the Lieb lattice has a magnetic ordering other than 
paramagnetism. To know which ordering it is, one needs 
to look at the results for q x and q y . 


The results for the magnetization amplitude to in the 
limit U —i 0 can be explained using first-order pertur¬ 
bation theory. Let us denote by Hq (the unperturbed 
Hamiltonian) the tight-binding terms of Eq. [ 6 ] that is, 
Hamiltonian Hmf with U = 0. Its eigenvalues are 


±2tyJ cos 2 + cos 2 


2 5 


±2t^ cos 2 (jf +q x ) +cos 2 + q y 'j, 


( 11 ) 


and two coincident flat bands at e = 0. Using the inter¬ 
action terms of Hamiltonian [ 6 ] as a perturbation yields, 
to first order, two key results. Firstly, the flat bands 
are split into two non-degenerate nearly flat bands. One 
of them is shifted to positive energy by an amount pro¬ 
portional to toU, while the other is shifted to negative 
energy, by the same amount, at each point of the Bril- 
louin zone. Secondly, the four non-flat bands are shifted 
by j-(riA + tlb)- These two conclusions allow us to pre¬ 
dict the behaviour of to near [7 = 0. For the following 
calculations, one must keep in mind that the diagonal 
term in Eq. [7] is also to be accounted for. 

We begin by filling up the lower bands (which cor¬ 
respond to the bands with minus signs in Eq. |11[ ), dis¬ 
tributing the particles among the sublattices according to 
Fig. 2b For n lower than 2/3, it is best to keep m = 0 


because, up to first order, the energy of the two lower- 
energy bands (associated with the Hamiltonian Hq) does 


C. Magnetic ordering 

Fig. [5] shows both the mean-field magnetic phase di¬ 
agram of the Lieb lattice (bottom plot) and that of the 
square lattice (top plot), for comparison. The phase di¬ 
agrams were computed using the independent variables 
n and U, in the range ( n,U ) G [0,2] x [0,20], and were 
obtained by joining the results for the three order param¬ 
eters: q x , q y and to. Near n = 0 and n = 2, the system 
is ferromagnetic for large U and paramagnetic for low U, 
like the square lattice, albeit with a wider ferromagnetic 
region. At intermediate U and n near 0.5 or 1.5, the 
system displays a (0, q^) spiral phase, characterized by 
q 2 « 7t/2. 

The spiral phase characterized by q = (n, n) only oc¬ 
curs at exactly n = 1, for any U, as in the square lattice. 
In the latter, this would be called antiferromagnetism. 
Nonetheless, in the Lieb lattice, a ( 7 r, n) phase should 
be identified with ferrimagnetic ordering uni] (see Fig. 
[3|. Indeed, the spin-spin correlation in a (7r, 7r) phase is 
ferromagnetic in each sublattice, but antiferromagnetic 
between different sublattices. The total spin per unit 
cell is finite, because to is finite at half-filling (see Fig. [4]) 
and the number of sites per unit cell is odd. 

When slightly doped away from half-filling (0.95 < 
n < 1.05), both q x and q y continuously deviate from 7 r 
and become a (< 71 , < 71 ) phase with <71 « 7 r. This area be- 
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FIG. 5. Mean-field magnetic phase diagrams of (a) the square lattice and (b) the Lieb lattice. In the case of the square lattice, 
the value of q varies continuously in the range [0, 7r], in each region labelled as such. In the case of the Lieb lattice, q\ ~ 7r and 
q 2 ~ 7t/ 2, and the transitions between regions labelled using qi or q 2 are discontinuous. 


comes narrower in the n direction as U grows larger. This 
phase can be interpreted as a ( 7 r — 6, n — S) phase with 
small S , that is, a local (looking at only a few unit cells) 
ferrimagnet with a slow modulation in the direction of 
spins along the lattice. At large U, when further doped, 
the system undergoes a first-order phase transition from 
q « (7r,7r) to 9 ~ (< 11 , 92 ) with q 2 ~ 7t/ 2, reflecting lo¬ 
cal antiferromagnetic correlations in the x direction, and 
sublattice-wise antiferromagnetism in the y direction. In 
other words, each sublattice is ferromagnetic in the x di¬ 
rection and antiferromagnetic in the y direction. If doped 
even further away from half-filling, two more first-order 
phase transitions occur: first to ( 91 , 0 ) and finally to ( 0 , 0 ) 
(ferromagnetism). At regular intervals in n (namely 0.11, 
0.22, 0.33 and their symmetric counterparts), we find fer¬ 
romagnetic dips into the paramagnetic region. These can 
most likely be explained using the symmetry of the lat¬ 


tice and higher-order corrections. 


IV. CONCLUSION 

In summary, we have computed and analysed the n- 
U mean-field magnetic phase diagram of the Lieb lat¬ 
tice, and compared it to that of the square lattice. Far 
from half-filling, the two phase diagrams display ferro¬ 
magnetism [9 = (0,0)] for high U and paramagnetism 
(to = 0) for low U, while at exactly half-filling (one elec¬ 
tron per lattice site) the ground state is a ( 7 r, 7 r) spiral 
phase for both lattices. 

Although the diagrams coincide at n = 1, it is close to 
that line that their most remarkable differences arise. In 
fact, at large U, as we move away from half-filling [(the 
(•7T, 7 r) phase)], the Lieb lattice undergoes three first-order 
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phase transitions (7r,7r) —»■ (7r, 7t/2) —> (7r,0) — > (0,0), 
unlike in the case of the square lattice, where the transi¬ 
tion from antiferromagnetism to ferromagnetism is con¬ 
tinuous in q. On the other hand, near the tight-binding 
limit and within the interval n £ [2/3,4/3], the magneti¬ 
zation of the Lieb lattice is finite and the ground state is 
characterised by spin spirals, contrasting with the para¬ 
magnetic ordering of the square lattice in this region of 
the diagram. 

Our numerical results are in agreement with a theo¬ 
rem by Lieb [2], which applies to Hubbard models which 
comprise sublattices with different number of sites (in 
our case, we have twice as many B/C sites as we have A 
sites). The theorem states that the ground state of such 
a system at half-filling is ferrimagnetic [IS]. According 
to our results, the ground state at half-filling is charac¬ 
terized by q = ( 7 r, n) and finite to, which translates into 
ferromagnetic sublattices and finite total spin on each 
unit cell (see Fig. |3|, which is qualitatively consistent 
with the theorem. According to this theorem, however, 
the total magnetization per unit cell in the Lieb lattice 
at half-filling should be 1 for any U. Our mean-field 
approach yields that value as U grows large but deviates 
from 1 at low U (see Fig. [4]). On the other hand, this the¬ 
orem is also applicable to a square lattice Hubbard model 
if one divides the lattice into two sublattices. This has 
been done before [26] with a square lattice divided into 
two sublattices, A and B, with the same number of sites 
each. In consonance with the aforementioned theorem by 
Lieb, this square lattice has zero total spin per unit cell 
at half-filling, for any U, regardless of to being zero or 
not. Therefore, it stands to reason to conjecture that the 


mean-field calculations performed for the square lattice 
also return wrong values for to at low U, even though the 
correct values cannot be deduced from Lieb’s theorem, as 
it only predicts the total spin per unit cell. 

The disparity between our mean-field results at half- 
filling and the prediction of Lieb’s theorem may be due 
to two important restrictions that we imposed in order 
to simplify our calculations. Firstly, we assumed that the 
occupation numbers for any U remain the same as in the 
tight-binding limit (U = 0), and secondly, we assumed 
that the magnetization is the same on every sublattice. 
If it turns out that these two assumptions are indeed the 
reason for the discrepancy, that is, if the Lieb’s theorem 
can be satisfied in a mean-field approach applied to this 
paper’s model, albeit with more relaxed constraints, such 
a result should be taken into account in any other mean- 
field study of interacting fermions in bipartite lattices or 
even more complex lattices whose unit cells contain more 
than two types of atoms. This is an open question that 
we intend to address in the future. 
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